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SYMBOLS 

a Amplitude of wall displacement 

A Fourier-Chebyshev amplitude of v; see (16) 

C Pressure coefficient 

P 

f Wall shape function 

G k See (20) 

H H(Y) = dY/dn ; see (17) 

k Wavenumber of wall displacement 

k Turbulent kinetic energy 

m X-wavenumber 

N Number of retained modes in or t, 
p Pressure 

Re Ux/v , Reynolds number 

t Time 

T Reynolds stress 

Chebyshev polynomial of degree n 
Eigenfunction; see (18) 
u Velocity in x-direction 

U Free-stream velocity in x-direction 

v Velocity in y-direction 

V k See (19) 

x Unmapped Cartesian coordinate 



X Mapped coordinate; see (14) or (15) 

y Unmapped Cartesian coordinate 

Y Mapped coordinate; see (14) or (15) 

a =1,2; coordinate label 

At Time step 

Ax Space discretization interval 

Au Velocity fluctuation 

5 Boundary layer thickness 

n Coordinate in conformal mapped system 

4> 0 Phase of wall displacement 


iii 


Eigenvalue; see (18) 
y T Eddy viscosity coefficient 
v Viscosity 

it Pi 

p Density 

p Parameter used in (17) 

a Parameter used in viscous stabilization; see (13) 
C Coordinate in conformally mapped system 

Superscripts 

Convective update 
Pressure update 
— Viscous update 

Viscous stabilization 

' Derivative with respect to £ + in ; 

Also, fluctuation from mean flow 



I. INTRODUCTION 


The study of mechanisms of possible drag 
reduction effects in flows over wavy walls as 
well as the study of mechanisms of generation 
of surface waves by wind requires understanding 
of the detailed flows over wavy surfaces. In 
this paper, we describe a computer code based 
on spectral methods to study two-dimensional 
incompressible flows in wavy geometries. 

Work on flows in wavy geometries began 
with the classical analysis of Kelvin and 
Helmholtz of linearized inviscid surface waves'*'. 

In the absence of mean- flow-shear effects, 
inviscid theory predicts that, for wall displace- 
ments of the form y = a cos (kx + $q) with 
ka <<1 (so that the problem may be linearized 
in terms of ka) , the pressure distribution is 
180° out of phase with the wall displacement 
and the pressure coefficient C = |p | /^pU^ 

has the value 2ka at the wall. Miles z extended 
the inviscid theory of account for shear in the 
mean flow, still restricting attention to very 
gentle waves. Miles' theory improves the pre- 
diction of the pressure coefficients. However, 
because of the inviscid nature of the analysis, 
the pressure phase shift at the wall is still 180° 
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Benjamin analyzed high Reynolds number laminar 
flows over wavy walls including the effects of 
shear still assuming ka <<1. Benjamin's theory 
typically predicts a pressure phase shift of 
210° relative to the displacement of a solid 
wavy boundary. The theory of Miles and Benjamin 
has also been extended to turbulent flows. Davis 
analyzes turbulent flow over wavy surfaces by 
postulating that Reynolds stresses are, to a 
first approximation, constant along lines of 
constant n, where (C,ri) is a first-order 
accurate orthogonal coordinate system with 
the wavy boundary at n= 0. Davis' analysis is, 
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therefore, also limited to small ka. Markatos 
has recently presented numerical results for 
transfer of heat, mass and momentum in flow 
over an evaporating wavy water surface. While 
Markatos' numerical results are at best first- 
order accurate in space, they go a long way 
towards providing understanding of the nature 
of flow over wavy walls. 

The work described in this paper extends 
previous work in several ways. The full 
time-dependent two-dimensional Navier-Stokes 
equations are solved using spectral methods 
to achieve high spatial accuracy and high-order 
time-splitting methods together with conformal 
mapping methods to allow simulation of flow 
over steep waves. Results for laminar flows at 
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Reynolds nuinbers up to Re = 4 * 10 are presented 

x 

in Sec. III. Turbulent flows over wavy walls are 
simulated by using a spatially varying eddy vis- 
cosity distribution. Some preliminary results 
and comparisons with experiment are given in 
Sec. IV. 


are 


II. METHODS OF SOLUTION 


The two-dimensional Navier-Stokes equations 


3 v 2 3 v „ 2 

a + Y v 2L = _ lE + Y _L_ T 

8=1 B9x B 9x ct 8=1 9x 8 aB 


3t 

with a =1,2, 


a 


2 3 v_ 

^ , 3x 
a=l a 


= 0 , 


(1) 


( 2 ) 


where the stress tensor may include both laminar 
(viscous) stresses, vVv,and turbulent Reynolds 
stresses (to be discussed briefly in Sec. IV below) 
Eqs. (1) - (2) are to be solved in the region 


0 1 X 1 < 27T ' f(x lf t) < x 2 < co (3) 

above the wall = f(x^,t). In the present 
paper, we restrict attention to surfaces of the 
form 

f(x^,t) = a cos kx (4) 

independent of t with periodic inflow-outflow 
boundary conditions in x^ , 


viz . 



v a (x i + 2ir ' x 2 ' 


t> = v a (x l' 


x 2 ' fc) 


(5) 


Extensions to study the effects of time-dependent 
geometries and inflow-outflow boundary conditions 
are presently being made. 

A conformal mapping technique^ is used to 
transform the region (4) into the region 

0 2tt t 0 < n < - 

If a £ 1, the conformal mapping coefficients are 
accurately generated in only order N £'og N opera- 
tions where N is the resolution along x-^ ( z, ) , 
so time dependent geometries can be efficiently 
handled. A further (non-con formal) stretching 
transformation of p is used to implement the 
spectral methods described below see Eqs. (14)- (15) 
Eqs. (1) - (2) are solved numerically by 

a fractional step procedure. The fractional steps 
are: (i) Convective update 

An intermediate field i s found so that 

v = v (t) + ^ AtN (t) - -kl (t-At)At ( 6 ) 
a a 2 a 2a 

where 

N = l v Q 3 -^— v 
a L 3 3v^ a 

Here a second-order accurate Adams-Bashforth scheme 
is applied. 
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(ii) Pressure update 


v 

a 


3g 3p _ 3 n 3p 

V a 3 x 85 3x 3n 

a a 


2 


z 

a=l 


A 


f 3 £ 3 v a 



+ 


A 


3 n 3 v, 

3x 3 n 
a 


0 


(7) 

(8) 


with the boundary condition 


2 


z 

a=l 



a 


0 at n =0. 


(9) 


Substituting (7) into (8) gives the Poisson 
equation 






3 n 3 

3x 3 n I 
a 1 


( 10 ) 


with 


F' 


2 


z 

a=l 



2 

z 


a=l '■ 


In ' 

3x 

a-' 



(ID 


where the conformal nature of the map from (x n , x 
to ( 5 , n ) is used. 



(iii) 


Viscous update 


v 

a 


v + 
a 


At 


T aB (t) 


with the boundary conditions 


(12) 


v =0 at n = 0. 
a 

This Euler step is only first-order accurate but 
will be improved upon later. The time step 
restrictions for numerical stability of (12) are 
particularly severe, especially in the case of 
turbulent flows, so the next fractional step 
is designed to stabilize the scheme. 

(iv) Viscous stabilization 


v - — 
a i T 


! 2 

1 min 


V 2 * * * * 7 v = 
? a 


v - 

a 


a 


If' 


2 

min 


v a (t > 


(13) 


with 

v = 0 at r| = 0 , 
a 

2 2 2 2 2 

where = 9 /9£ + 9 /9n is the Laplacian 

with respect to £,ri . Here a is a parameter 
chosen to achieve optimal stabilization and accuracy, 

(v) Local extrapolation 

For time dependent problems requiring 
accurate calculation of transient effects, a local 
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extrapolation procedure xs used to achieve global 
second-order accuracy in time. Further details and 
caveats on the splitting procedure described above 
are given in Ref. 7. 



In order to implement this algorithm in a 

g 

spectral method , the semi-infinite region 0 £ 5 _< 2tt, 
0 _< p < <=° is mapped into the finite rectangular domain 

0 < X < 2tt, -1 <_ Y < 1 
by 


X = S , 



(14) 


or 

X = ? , Y = 1 - 2e" n/L (15) 

for some suitable mapping scale L. 

Next, all dependent variables are expanded 
in Fourier-Chebyshev series of the form 

N . 

v (X,Y,t) = l 7 A (n,m,t)e imX T (Y),(16) 
n=0 |mf<M 

where (Y) = cos (n cos ^ X) i s the Chebyshev 

polynomial of degree n. 

Technical details regarding spectral methods 
and their implementation are given in Ref. 8 
and references given therein. For illustration 
purposes, the method of solution of the Poisson 
equation (10) or (13) will now be described. After 
the transformation (14) or (15), the equation to 
be solved is of the form 

2 

pu + + H (Y) ^|h(Y)|£ = g (X, Y) . (17) 
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Introducing the eigenvalue problem 


H <Y) 5 $ H(Y) 55 - = X k u k (18) 

subject to the appropriate homogeneous boundary 
conditions, say u^ = 0, it follows that 


U (X , Y) = 

l V X) V Y > 

k 

(19) 

g(x,Y) = 

l G k (x)u k (Y > 

k 

(20) 


where V, satisfies 
k 


P V k + ~T 


+ A k V k = G k * 


Then (21) is solved by Fourier transformation in X, 
so the solution u to (10) or (13) can be found 
by further matrix operations. 

In the next Section, results are reported 
from this code using up to 64 Fourier modes and 
Chebyshev polynomials in each of the X and Y 
directions. Typically, the fractional step method 
with local extrapolation and tensor stress evalua- 
tion requires about 6 ms per retained mode per 
time step on a CDC 6600 computer. Time steps are 
restricted for stability to be less than the 
convective stability limit Ax/U . 
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Ill 


LAMINAR FLOW RESULTS 


Results from three laminar flow runs will 
be described here. The first case is flow over a 
very gentle wave a cos kx with k = 4-rrcm ^ , 

5 

ka = tt/ 1000 at a Reynolds number Re = 4 x 10 

X j - . 

at x = 29 cm, so the boundary thickness <5 = 6.5/ 

= 0.3 cm. A plot of the mapped grid used in the °° 
calculations of this flow is given in Fig. 1, 
where only the portion of the grid near the wall 
is plotted. Since the x and y scales are different, 
the conformal nature of the map is obscured in 
the plot. 

In Fig. 2 we plot contours of the pressure 
at t = 90 (in nondimensional units in which U = 1) 
after the solution has converged to a steady state. 
The calculation was performed using 32 Fourier 
modes in X and 33 Chebyshev polynomials in Y. 

The contouring is done in the unmapped x^ = x 
and = y coordinates. The phase of the pressure 
distribution lags 214° ±5° behind the surface 


displacement a cos x and the pressure coefficient 

I I /I tt2 • __ 


C = 
P 


,Ap u; 


C ^ 1.76 x 10 
P 


at the wall. As mentioned above, the Kelvin-Helmholtz 

. . . . —3 

inviscid analysis predicts Cp=2ka - 6.3 x 10 

On the other hand, Benjamin's theory 3 predicts a phase 

shift of about 217° with in good agreement with 

our numerical simulations. 
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In Figs. 3 and 4 we plot contours of the velocity 
in the direction in the unmapped x^ = x 

and x -2 = y coordinates. Fig. 4 is a blown-up 
portion of Fig. 3 that highlights the structure 
of the velocity near the wall. Note that the 

scale of the boundary layer is such that u = 0.9611^ 
at ky - 2.5. In Fig. 5 we plot contours of u^ 

in the x-y plane blown-up to highlight the wall 
region. In Figs. 6 and 7, we plot contours of the 
spanwise vorticity at t = 90. Again, Fig. 7 is 
a blown-up version of Fig. 6. 

The second run is similar to the first except 
that the amplitude of the surface displacement is 
increased to 

ka = 5.5 x 10 -2 


In Fig. 8, a contour plot is given of the pressure 
distribution at a nondimensional time t = 90. 

The phase of the pressure distribution now lags 
about 212° behind the surface while the pressure 
coefficient is increased to about 

C - 3. 3 x 10 -3 . 

P 

Benjamin's theory still predicts a phase shift of 
about 217° with in good agreement with the value 
obtained by numerical simulation. In Figs. 9-11, 
we give contour plots in the wall region for the 
velocity components v ^/ v 2 and spanwise vorticity, 
respectively, in the second run. 


The third run involves the simulation of flow 

g 

studied experimentally by Kachanov et al . Here 
k = 2.5 cm ^ while 

ka = (5.3 ± 0.1) x 10 -2 . 


The calculation was performed at a Reynolds number 

5 

Re = l.l x 10 to match the experiment, while 

U ro = 530 cm/s so the boundary layer thickness 

/ ' 

77 — at x = 3.06 cm is 6 = 0.6 cm. 

Because the wavelength of the wall displacement 
is much longer in this run than in the second run, 
it turns out that the accurate simulation of this 
case is somewhat more difficult. Pressure pertuba- 
tions due to the wall die out over a distance of 
order the wall wave-length which is much larger 
than 6 so the resolution requirements in Y 
are more severe than for the first two cases. The 
results reported below were obtained using 32 
Fourier modes in X and 65 Chebyshev polynomials 
in Y. 

In Fig. 12, we plot contours of the pressure 

field for the third run after a steady state has 

been achieved. In Figs. 13 and 14, contour plots 

of the spanwise vorticity and velocity v^ , 

respectively, near the wall are given. Kachanov 
9 

et al measured the distribution of x-velocity 
v^ = u fluctuations at a height of about 6 above 
the wall. They found a fluctuation amplitude Au of 


Au 


e xpt 


U 


OO 


0.005 ± 0.002 
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with a phase shift of 40° ± 5° relative to the wall. 
The computational results give 
Au 

— S2ER = 0.004 ± 0.0002 

^OO Q ^ 

with a phase shift of 22 ± 5 . The magnitude of the 
fluctuations agree well between laboratory and 
computer experiment. The discrepancy in the phase 
shifts may be due to the sensitivity of this quantity 
to measurement error or due to the precise manner 
of measurement (in the computer experiment, the 
fluctuations of u were measured at constant Y while 
in the laboratory experiment the measurements 
seem to have been conducted at constant y) . 


IV. TURBULENT FLOW RESULTS 

In this Section, we report a simulation of 
Kendall's experiment" 1 of turbulent flow over a 
wavy wall. An eddy viscosity model is used to 
simulate the effects of the turbulence. The turbulent 
Reynolds stresses are modelled as 


.2 

pv^ = 

2 “ 

; 2y T 3v 1 /3x 1 - j pk 

(22) 

pv i v 2 

= y T (3v 1 /3x 2 + 3 v 2 /3x 1 ) 

(23) 

. 2 

P v 2 

= 2y T 3 v 2 /3x 2 - j pk 

(24) 
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where primes indicate turbulent fluctuations, overbars 
are averages, y T is a (spatially varying) eddy 
viscosity coefficient and k is the turbulent kinetic 
energy. 

The eddy coefficient y is modelled by use 

1 11 

of a turbulent boundary layer code that includes 

effects of wall curvature and pressure gradients 

12 

within a one-equation eddy- viscosity model. The 
actual pressure measurements of the Kendall experiment 
were used in this boundary layer code. The turbulent 
kinetic energy k is modelled on the basis of flat 
plate data according to 


- vjv' / k « 0.15 (25) 

A contour plot of the resulting eddy viscosity 
coefficient y ^ is given in Fig. 15. 

j. 

Kendall's experiment is simulated using a 

wall displacement a cos (kx) with a = 0.32 cm and 

k = 0.6 2 cm . The mean flow velocity in the free 

stream is U = 500 cm/s while the molecular viscosity 
2 

is 0.14 cm' /s. The calculations are performed using 
32 Fourier modes in X and 33 Chebyshev polynomials 
in Y. The flow field is initialized using the result 
obtained from the boundary layer code'*"'*' corrected 
to be incompressible. 

In Fig. 16, we plot contours of the velocity 
field v-^ in the x^-direction at t = 2s. In Fig. 17, 
spanwise vorticity contours at t = 2s are shown. 

Finally, in Figs. 18 and 19, contour plots of v^ and 
V 2 , respectively, at t=4s are given. It can be observed 
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from Fig. 18 that a small separation bubble is developing 
in the flow at t=4s. Separation seems to begin near 
t=3.2s. While no separation was observed in Kendall's 
experiment, it is known that the laboratory flow is 
close to separation. There are several possible reasons 
for the appearance of a separation bubble in the 
computation. First, the eddy viscosity distribution y T 
is fixed for all time by the turbulent boundary layer 
calculation. Since the mean-flow velocity has evolved 
significantly from the boundary layer approximation, 
it seems reasonable to expect that u T should change 
also. Second, the modelling of normal Reynolds stresses 
in (22)- (24) may be suspect. A possible remedy for 
both these problems would be integrating a multi- 
equation turbulence model at the same time as the mean 
flow equations. Further work is now planned in this 
direction. Finally, while the numerical resolution of 
the flow seems adequate in the wall region, there does 
seem to be need for more resolution at large distances 
from the wall. For this reason, higher resolution 
calculations with modified mappings in the y direction 
are now underway. 


V. CONCLUSION 

A computer code to solve the full two-dimensional 
Navier-Stokes equations in an arbitrary wavy geometry has 
been developed. The code has been run on both laminar 
and turbulent flows over wavy walls. The results have 
been compared with both available theory and experiment. 

For laminar flow, good agreement is achieved. For turbulent 
flow, the agreement is less satisfactory and may be explained 
by the simplified turbulence modelling used in these 
calculations . 
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Figure 1. A plot of the coordinate system near the wall 

for Run 1. Here the wall displacement is y = a cos kx 

-4 

with k = 10 tt and a=10 
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Figure 4. Same as Figure 3 except that the region 
near the wall is expanded. 






Figure 7 . Same as Figure 6 except that the region near 


the wall is expanded. 
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Figure 18. Contour plot of 

of Kendall's experiment. 


at t=4s in the simulation 
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